Disks around Young Planetary-mass Objects: Ultradeep Spitzer Imaging of NGC 1333

NGC 1333 was observed as part of the YSOVAR program during Spitzer's Warm Mission, under program id 61026. The field was observed in IRAC channels 1 (3.6 μm) and 2 (4.5 μm), in the following called IRAC1 and IRAC2. For full information on YSOVAR, see Rebull et al. (2014). In short, YSOVAR was a large-scale monitoring program with Spitzer to study variability in young stars, focused in particular on inner disk structure, accretion changes, eclipsing binaries, and rotation periods. The program surveyed 12 star-forming regions. The number of epochs varies from region to region, from a minimum of about a dozen to several thousand. NGC 1333 was one of the target regions for YSOVAR, with a survey area of 2 × 2 IRAC field of views, where a $10^{\prime} \times 10^{\prime}$ field was covered by both channels. For more information on the variability analysis of the NGC 1333 data set, see Rebull et al. (2015). The fundamental reference for the IRAC instrument is Fazio et al. (2004).

We downloaded 70 images from the YSOVAR data set on NGC 1333, for each of the two channels, from the Spitzer Heritage Archive. The images for a given channel all have very similar pointings, with minimal field rotation, and similar depth. However, due to the layout of the Spitzer focal plane, the images for channel 1 are offset compared to channel 2. Stacking was performed using the Python reproject package (Robitaille et al. 2020), with the function reproject_and_coadd. We stacked in two iterations, for each filter separately: first we created seven stacks from 10 images each, then we stacked the resulting seven stacks to a final deep image for each band. In Figure 1 we show a comparison for a part of the stacked image versus the same part of a single image.

Zoom In Zoom Out Reset image size

We created a catalog of brown dwarfs in NGC 1333 based on the recent census by Luhman et al. (2016), which builds on several previous surveys and includes a comprehensive literature review. Their published list comprises 203 new members, the clear majority of which are in the area covered by our deep Spitzer images. Most members have optical or near-infrared spectroscopic confirmation, either by the authors of that census, or from the literature. In addition, the authors use their own astrometry, as well as indicators of youth (e.g., X-ray emission, infrared excess) to identify cluster members. The Luhman census is "nearly complete" down to K magnitudes of 16.2 and an A_J of 3.

From that list, we selected the 65 for which the spectral type adopted by Luhman et al. (2016) is M6 or later. A spectral type of M6 is commonly used as the boundary between stars and brown dwarfs in star-forming regions. It translates to a temperature of ∼3000 K (Mužić et al. 2014), corresponding to an object with a mass at or slightly above the substellar limit of 0.08 M_⊙ for an age of 1–5 Myr, according to evolutionary tracks (Baraffe et al. 2015). That means our sample should safely encompass all known brown dwarfs in the NGC 1333 census, but may include some very-low-mass stars as well. The sample also includes 19 objects with a spectral type of M9 or later, i.e., with estimated masses near or below the deuterium-burning limit. Those PMOs are the main focus of this paper.

Of the 65 likely brown dwarfs, 49—including 14 PMOs—are covered by at least one of our deep Spitzer/IRAC images. After the Luhman et al. (2016) census, two more papers confirmed new substellar members for NGC 1333: Esplin & Luhman (2017) and Allers & Liu (2020). The former includes two more objects with M9 or later spectral type, one of which is covered by our images, and was originally reported by Oasa et al. (2008). We added it to our sample (with spectral type M9–L4, we adopt L1.5 here). Thus our total sample of PMOs is comprised of 15 objects. Due to the spatial offsets between the images in the two Spitzer/IRAC channels, a few objects in our sample are only covered by one channel. A dozen of the sources are in areas with a highly variable background or very close to the edge. Of the 65 likely brown dwarfs, 34 are in the JKS catalog published by Scholz et al. (2012), which combines the deep J- and K-band photometry from Scholz et al. (2009) with the list of young stellar objects created from the the single-epoch Spitzer/IRAC images taken by the C2D project (Evans et al. 2009). For the objects covered by our deep Spitzer images, 23 (IRAC1) and 22 (IRAC2) are in the JKS catalog.

We characterize our sample in Figure 2. In the left panel, we show spectral type versus K-band magnitude for the whole brown dwarf catalog, marking those covered by our deep images, and indicating approximate mass limits. As can be appreciated from this figure, our sample covers the entire substellar range, from the limit between stars and brown dwarfs to masses well below the deuterium-burning limit. In particular, it includes all objects with spectral types L0 or later found in this cluster thus far. In the right panel, we plot the spatial distributions of the brown dwarfs and PMOs (marking those with disks; see Section 3). The figure demonstrates that our images cover the embedded cluster as well as the surrounding areas.

Zoom In Zoom Out Reset image size

We carried out aperture photometry for the objects in our catalog using tools from photutils (Bradley et al. 2021). We chose a constant aperture with a radius of 5 pixels. We note that the native pixel size of IRAC is 12, but the mosaics used here have 06 per pixel. The FWHM of IRAC during the Warm Mission is 20, according to the IRAC Instrument Handbook. An aperture radius of 5 pixels corresponds to an aperture diameter of 6'', or approximately 3 times the FWHM, and should capture >90% of the flux. For local background estimates, we measured the median flux in an annulus with an inner and outer radius of 7 and 10 pixels, respectively, centered on the object coordinates.

The magnitudes were calibrated by direct comparison with the JKS catalog mentioned above. The offset between the JKS IRAC magnitudes, which originally come from the C2D source list (Evans et al. 2009), and our instrumental measurements is constant for a wide range of magnitudes, with standard errors of 0.03 and 0.04 mag (depending on the channel). Thus, simply adding an offset shifts our instrumental magnitudes into the standard system and makes them comparable with the literature.

Photometric errors were estimated by adding in quadrature the Poissonian errors in the source flux, the standard error for the background times the number of pixels in the aperture, and the error in the calibration. The median error, which is dominated by the calibration uncertainty, is 0.04 mag for both channels.

All sources were checked in the stacked images. In particular, we verified that the background determined by the photometry routine is adequate and that the aperture is well centered on the source. As indicated earlier, some objects are located in areas of highly variable background, which may affect the photometry.

Finally, we calculate the K – IRAC1 and K – IRAC2 colors using the K-band values listed in Luhman et al. (2016). We dereddened these colors to zero extinction using the published A_J extinction values for the sample, again from Luhman et al. (2016). Those values have been determined by comparing the optical and near-infrared colors to the expected photospheric values. For a subset of our sample, Scholz et al. (2012) determined extinctions, using photometry and spectroscopy. Comparing with the Luhman values, there is good agreement, with a typical uncertainty <0.5 mag in A_J .

For the extinction correction, we convert from A_J to A_K using a standard extinction law (A_K = 0.382 A_J ; Mathis 1990). For the IRAC wavelengths, we adopt A_IRAC1 = 0.6 A_K and A_IRAC2 = 0.5 A_K . These extinction relations are in line with recently published work for star-forming regions; see Chapman et al. (2009). With very few exceptions, the brown dwarfs in our sample have an A_J below 5. Thus, the exact choice of the extinction law within plausible values does not change the outcomes in any significant way. The dereddened colors versus spectral types are plotted in Figure 3, as a way to identify objects with infrared excess emission and thus disks. These figures are the foundation for the discussion in Section 3.

Zoom In Zoom Out Reset image size

For completeness, we also did full-field photometry using the SExtractor (Bertin & Arnouts 1996) for the stacked image and an individual image in both channels. Comparing the peak of the histogram for the magnitudes, we find that the stacked images are deeper by 2.4 and 2.6 mag, for IRAC1 and IRAC2, respectively. This result confirms the face-value expectation—stacking 70 images of similar quality should result in $\sqrt{70}=8.3$ better signal-to-noise ratio, or 2.3 mag greater depth.

To determine which objects have excess emission in their infrared colors we first need to establish the photospheric levels in the same colors, as a function of spectral type. We prepared a sample of stars and brown dwarfs with spectral types from M0 to L1, all without disks (and thus without infrared excess), that are members of nearby star-forming regions (<10 Myr). In this exercise, we are building on the work in Almendros-Abad et al. (2022). The K – IRAC1 and K – IRAC2 colors for this sample show a gradual increase from values of 0.1 for early M dwarfs to around 1.0 for early L dwarfs. We fit these colors for spectral types M6 or later with a linear slope, which describes the trend seen in the data well:

$$\begin{eqnarray}&&K-\mathrm{IRAC}1=0.1007\times \mathrm{SpT}-0.2847\end{eqnarray} \tag{ 1 }$$

,

$$\begin{eqnarray}&&K-\mathrm{IRAC}2=0.1179\times \mathrm{SpT}-0.3281\end{eqnarray} \tag{ 2 }$$

.Here, the spectral type (SpT) is 6.0 for M6 and 11.0 for L1. We tested these relations by doing the same procedure on only a subset of the diskless objects, and do not find a significant change. In particular, we limited the fit for objects with ages <5 Myr and the fit remains very similar. For reference, the given relations are consistent with the ones published by Luhman (2022) for ages <20 Myr.

Our photospheric colors are plotted in Figure 3 as solid lines. We also show in the figure the typical scatter of diskless objects around this fit, as the colored area. Objects without infrared excess in NGC 1333 are expected to be found in that area around the solid line in Figure 3.

In Table 1 we list the newly derived infrared photometry and colors (before extinction correction) for the 15 PMOs (spectral types M9 or later): 14 have measurements in IRAC1, 13 in IRAC2. Three objects with reported magnitudes are in areas of highly variable background, and after careful examination we deem those measurements to be unreliable. As can be seen in Figure 3, eight objects with good measurements have a dereddened K – IRAC1 color consistent with photospheres, six do in K – IRAC2. Altogether seven sources out of 12 with robust measurements in at least one band show no evidence for infrared excess.

The remaining five objects are found to have infrared excess and are thus likely to be PMOs with disks. Four of these are at spectral type M9, three of them with excess in both bands (one with marginal excess in IRAC1). The fifth, with spectral type L0, has no excess in IRAC1, but a clear excess in IRAC2. Thus, the disk fraction from our sample is 5/12 or 42%. Due to the small sample size, the plausible range for the disk fraction is 26–60% (assuming 1σ binomial confidence intervals). Visual examination of the three objects in areas with highly variable background shows that one (SONYC-NGC1333-9) is bright in both IRAC bands, and thus also a good candidate for harboring a disk. The other two are faint and more likely to lack excess. Including those three would bring the disk fraction to 6/15, or 40%. We note that with one exception all disks in the PMO regime are found for brown dwarfs classified as M9, at the adopted upper limit for our PMO sample. For L0 or later objects, the disk fraction is only 1/5 or 20%, 1/6 if we add the one with background-affected measurements.

Our analysis includes the lowest-mass objects identified thus far in this cluster, in particular, the six objects with spectral type L0 or later (corresponding to masses of ∼0.01 M_⊙ or lower). We determine the infrared excess for five out of six objects, with the one exception mentioned above. It is likely that the current optical/near-infrared surveys are not complete for those spectral types. Therefore, the disk fractions determined here may be affected by incompleteness at the lowest masses. If disk-bearing objects are on average deeper embedded than those without, thus having more extinction and are fainter, it is conceivable that we are underestimating the disk fraction for M < 0.01 M_⊙. This, in addition to small number statistics, may contribute to the the low fraction of disks in the L0 or later domain. Deeper observations will be needed to clarify this.

In Figure 2, right panel, we show the spatial distribution of the PMOs and PMOs with disks relative to the cluster extent and the overall brown dwarf sample. There is no obvious spatial bias in those samples, they are all concentrated on the cluster itself. One of the disk-bearing PMOs is located just outside the nominal cluster extent, but the same is true for several PMOs without disks.

Among the brown dwarfs with valid measurements and spectral types M6 to <M9 in our images, 8/27 show clear infrared excess in IRAC1, 15/23 in IRAC2. This corresponds to a disk fraction of 30% ± 10% for IRAC1 and 65% ± 12% for IRAC2. Thus, based on our measurements, the disk fraction among PMOs is consistent within the statistical uncertainties with the value for more massive brown dwarfs. We only quote these numbers for completeness: in contrast to the faint PMOs, our deep images do not provide a significant benefit for the brighter M6 to <M9 objects.

In Table 1 we also included the previous classification of the infrared excess based on Luhman et al. (2016), in the column labeled "IR". This previous paper uses mid-infrared excess measured from Spitzer images to identify circumstellar disks. Of our list of 15, 12 have that information listed in Luhman et al. (2016), and for all 12 our classification is in agreement with the one in the literature. We also classify three objects for the first time here. One of them likely has a disk: at spectral type L0, it is currently the latest type object with a disk in NGC 1333. We also classify the three objects with spectral types later than L0 as diskless, two of them for the first time.

Luhman et al. (2016) also derived disk fractions for subsamples divided by spectral type. They report 11/21 or 52% for spectral types later than M8 (a slightly different subset than our sample of PMOs), 20/33 or 61% for M6–M8, 28/48 or 58% for M3.75–M5.75 and 27/39 or 69% for K6 to M3.5. Within the statistical error bars, these values are all consistent with each other, and with our own measurements of the disk fractions, given above. Our measurement is also consistent with the disk fraction of 66% ± 8% for objects with spectral type of M5–M9 in NGC 1333, reported by Scholz et al. (2012). For completeness, for a sample of 79 Gaia-selected low-mass stars in NGC 1333, Pavlidou (2022) measure a disk fraction of 67% ± 13% (for the mass range 0.1 to 1.0 M_⊙, approximately).

We note here that the "disk fraction" in the quoted papers and in our study is defined as the number of Class II objects, divided over the sum of Class II and Class III. These numbers are all based on samples from optical/near-infrared surveys, and should not be compared directly with the higher disk fractions derived from mid-infrared surveys (Jørgensen et al. 2006; Gutermuth et al. 2008). One further caveat: disk surveys based on photometry at 3–8 μm can by their nature not find disks with large inner holes where the disk only causes excess emission at longer wavelengths. All studies mentioned above share this bias.

In Figure 4 we provide a summary of the currently available disk fraction measurements for low-mass stars, brown dwarfs, and PMOs in NGC 1333, including only objects with known spectral type. As can be appreciated from this figure, there is no significant trend in the disk fraction from late K to early L spectral type. In that regard, our new photometry confirms results from previous studies of the disk population for this cluster. Finding a disk fraction largely independent of mass for low-mass stars and brown dwarfs is also in line with what has been found for other star-forming regions with ages between 1 and 5 Myr, including IC348 (Luhman et al. 2016), σ Orionis (Scholz & Jayawardhana 2008), and Chamaeleon-I (Luhman et al. 2008). As a sidenote, the brown dwarf disk fraction derived in the current paper from IRAC1 excess (30%) seems to be an underestimate compared with the literature and should be treated with caution.

Zoom In Zoom Out Reset image size

Figure 4 also illustrates that there may be a decline in the disk fraction when comparing L-type PMOs to K–M stars and brown dwarfs. As noted above, among the six cluster members in NGC 1333 with spectral type L0 or later only one seems to have a disk, based on our measurements, translating to a disk fraction of 20%. Given the small sample size and possible bias, as mentioned above, it is premature to conclude that this is indeed a drop-off in the disk fraction. If confirmed, it would constitute a very curious result, with possible implications for the formation process. Finding more L-type objects in this cluster would immediately give us a better handle on this issue.

We note that for somewhat older star-forming regions, like Upper Scorpius and the TW Hydrae Association, with ages ∼8–10 Myr, there is some evidence in the literature for an increased disk fraction for brown dwarfs compared to low-mass stars (Riaz & Gizis 2008; Luhman & Mamajek 2012; Cook et al. 2017). For those regions, PMOs do not seem to stand out from more massive brown dwarfs in terms of their disk fraction—but the small sample size hampers a more definitive assessment of the longevity of disks around PMOs.

Generally speaking, the evidence is strong that free-floating PMOs can harbour disks, and therefore, have the potential to form their own (miniature) planetary systems. In terms of the overall mass and scale, these systems would be more comparable to Jupiter's Galilean moons than to the solar system.

Our work in NGC 1333 reported in this paper establishes that at least five PMOs in NGC 1333 show evidence for the presence of a circum-substellar disk. The lowest mass among these is an object at spectral type L0 (not previously identified as disk bearing), corresponding to a temperature of 2200 K, which for an age of 1–2 Myr would result in an estimated mass of ∼0.01 M_⊙ (Baraffe et al. 2015). The three confirmed sources with later spectral types, and thus presumably lower masses, do not have detectable infrared excess and are classified as diskless.

Similar studies have been carried out in other star-forming regions. Taken together, there is now a small sample of young PMOs with securely identified disks. This includes OTS44 (Joergens et al. 2013), the only one with an ALMA detection at far-infrared/submillimeter wavelengths (Bayo et al. 2017). For young PMOs, definitively detecting an infrared excess due to a disk is challenging work, for two reasons. One, as can be appreciated from Figure 3, the photospheric near-/mid-infrared colors increase steadily for young late M and L dwarfs. With the typical uncertainties in the spectral typing of at least ±1 subtype, the photospheric level for a given source is uncertain. Second, these are faint sources, and thus, the near-/mid-infrared magnitudes come with considerable errors. Both difficulties worsen toward later spectral types and lower masses. They can be overcome with deeper mid-infrared images (the motivation for this study), improved spectroscopic classification, or photometry at wavelengths where the photospheric flux is negligible compared to the disk excess (10 μm or beyond).

With all that in mind, we summarize the results for regions with sufficient survey depth.

In σ Orionis, two planetary-mass brown dwarfs (SOri 60 and SOri 71) have consistently been found to host disks at wavelengths up to 8 μm (Luhman et al. 2008; Scholz & Jayawardhana 2008). These two have spectral types of L2 and L0, respectively. This cluster is slightly older than NGC 1333; for its age of 3–5 Myr, these spectral types would correspond to ∼0.01 M_⊙ or above. There are some PMOs with later spectral types and claimed excess emission (e.g., SOri 65, 66, and 70), but the mid-infrared data are inconclusive (Luhman et al. 2008; Scholz & Jayawardhana 2008). In the case of SOri 70, it is also not clear yet if it is in fact a young member of the cluster and not a foreground brown dwarf (Zapatero Osorio et al. 2017).

In Chamaeleon-I, OTS44 (M9.5) is a secure detection. In addition, Luhman et al. (2008) found Cha 1107–7626, an L0 object with a disk. Esplin et al. (2017) report three more objects with M9–L2/L3 spectral types that may have excess emission, but caution against interpreting it as evidence for a disk. Given the uncertainty in their spectral type, more work needs to be done before assigning masses or infrared excess to those.

In Taurus, the lowest-mass objects with safely detected disks have late M spectral types, according to the recent survey by Esplin & Luhman (2019). Some L0–L1 sources also may have infrared excess. In addition, Esplin & Luhman (2019) find one very-low-mass source with an L3 spectral type, which would have a mass well below 0.01 M_⊙, with possible excess emission at 4.5 μm. Again, this is not a robust detection yet.

ρ Ophiuchi has spectroscopically confirmed brown dwarfs with disks at spectral type M9–L1 (Testi et al. 2002; Jayawardhana & Ivanov 2006; Alves de Oliveira et al. 2012), which again would correspond to a mass of ∼0.01 M_⊙. CFHTWIR-Oph-33 has a spectral type of L3–L4 and may have excess emission, but the mid-infrared data are inconclusive in that regard (Alves de Oliveira et al. 2012; Esplin & Luhman 2020). Some other very-low-mass objects discussed in these papers may have infrared excess, but also very uncertain spectral type (e.g., CFHTWIR-Oph-58, with a range from M8.5 to L3; see Almendros-Abad et al. 2022).

Finally, the older Upper Scorpius association has a large population of brown dwarfs. Two with a spectral type of L3.5 are reported to have accretion and infrared excess (Lodieu et al. 2018). For ages of 5–10 Myr, this would correspond to 0.01 M_⊙ or above. To our knowledge, no lower-mass objects with disks have been found here.

In summary, at this stage, the low-mass limit for objects with a safe disk detection, across all star-forming regions, is ∼0.01 M_⊙ (or ∼10 M_Jup). The surveys have found perhaps two dozen objects in total that may have masses below that limit (depending on how the mass is estimated) in all of the very young regions mentioned above. For the majority of them we can already firmly rule out the presence of disks—the lowest-mass objects in NGC 1333 discussed in this paper belong in that category. For about ten of these PMOs with putative masses <0.01 M_⊙, the currently available data are inconclusive and insufficient for proper characterization. Therefore, it is certainly too early to rule out the presence of disks around objects with masses below 0.01 M_⊙. Having said that, given that the disk fraction among brown dwarfs is 30%–60% for the young regions mentioned above, the lack of confirmed disks below the 0.01 M_⊙ limit is perhaps already worth noting, especially if taken together with our new measurements for NGC 1333, as illustrated in Figure 4. If the disk fraction were similar for objects with masses <0.01 M_⊙, essentially all objects with tentative detections of infrared excess should have disks. Perhaps we are beginning to see a decline of disk fractions in the planetary-mass domain.

Sensitive mid-infrared observations of known PMOs, and deep surveys for new objects in this mass domain are needed to verify this suspicion. JWST is primed to help with both of these tasks. In NGC 1333 specifically, a Guaranteed Time program with NIRISS/WFSS is dedicated to a deep spectroscopic survey, and designed to find PMOs in this cluster down to masses of 1 M_Jup (program ID 1202; see Willott et al. 2022). Also, the unprecedented sensitivity of the MIRI instrument will enable us to probe directly for infrared excess in the lowest-mass free-floating objects known, at wavelengths exceeding those of Spitzer/IRAC, allowing for unambiguous disk detections.

The presence of disks is discussed in the literature as one possible signature to distinguish between a star-like and a planet-like formation scenario (Luhman 2012; Testi et al. 2016). In particular, if an object gets ejected early in the evolution from a forming planetary system, either by planet–planet scattering or encounters with other stars, it can be expected that its disk, if present, will be affected. Simulations on this issue, however, remain sparse or nonexistent. For the case of planet–planet-scattering, Bowler et al. (2011) demonstrate that the disruption of a circumplanetary disk should be common. In this context, disks detected with ALMA around young giant planets on wide orbits are an interesting comparison. These planetary-mass companions may have very compact disks (Wu et al. 2017), but it is unclear whether or not their masses are lower than expected from the standard relation between disk mass and stellar mass (Wu et al. 2020). ⁷

Given the discussed arguments, it is plausible to assume that free-floating planets that have been ejected from a young planetary system would not host massive, long-lasting disks. If future observations confirm the absence of disks for objects with masses below 0.01 M_⊙, it may indicate that this is the low-mass limit for objects to form like stars, physically set by the opacity limit for fragmentation.